function [qvalue,pvalue] = cusum_random(epsilon, direction)
    n = length(epsilon);
    X = zeros(1, n);
    S = zeros(1, n);
    epsilonList = epsilon;
    for i = 1:n
        X(i) = 2 * epsilonList(i) - 1;
    end

    if direction == 0 % forward cumsum
        S(1) = X(1);
        SMax = abs(S(1));
        for i = 2:n
            S(i) = S(i-1) + X(i);
            if abs(S(i)) > SMax
                SMax = abs(S(i));
            end
        end
    else % backward cumsum
        S(1) = X(n);
        SMax = abs(S(n));
        for i = 2:n
            S(i) = S(i-1) + X(n-i+1);
            if abs(S(i)) > SMax
                SMax = abs(S(i));
            end
        end
    end

    sum1 = 0.0;
    for i = floor((1 - n/SMax)/4) : floor(((n/SMax) - 1)/4)
        sum1 = sum1 + normcdf((4*i + 1) * SMax / sqrt(n)) - normcdf((4*i - 1) * SMax / sqrt(n));
    end

    sum2 = 0.0;
    for i = floor((-3 - n/SMax)/4) : floor(((n/SMax) - 1)/4)
        sum2 = sum2 + normcdf((4*i + 3) * SMax / sqrt(n)) - normcdf((4*i + 1) * SMax / sqrt(n));
    end

    pvalue = 1.0 - sum1 + sum2;
    qvalue = pvalue;
    
    %result = struct('p', pvalue, 'q', qvalue);
end